% task6_2_3_2.m

clear;
ConstantDefinition();
load("..\data\ConstantDefinition.mat", "FRE_F", "SampRate");
load("..\data\task6_2_2_4.mat", "sortP", "sortf", "Fs", "T", "fmt");

% 尝试自动化，一次失败的尝试
% funds = cell(size(sortf,2),1);  % to save fundamental results
% fPs = cell(size(sortf,2),1);
% Pth2 = 1e-5;
% 
% for i = 1:size(sortf,2)
%     % fundamental = [];
%     col = sortP(:,i);
%     if max(col)>Pth2
%         f = sortf(1:5,i);       % 5 most powerful components
%         fp = sortP(1:5,i);
%     else
%         f = [];
%         fp = [];
%     end
%     if isempty(f)
%         f = sortf((col>Pth2),i);
%         fp = sortP((col>Pth2),i);
%         if isempty(f)
%             f = [];
%             fp = [];
%         end
%     end
% 
%     [fundamental, funP] = ffilter(f, fp);   % detect fundamental component
%     % fundamental may equals [1312,656], as f is power-sorted
%     % so we must sort fundamental and ffilter again
%     [sf, si] = sort(fundamental,'ascend');
%     sp = funP(si);
%     [fundamental, funP] = ffilter(sf, sp);
% 
%     funds{i} = fundamental;
%     fPs{i} = funP;
% end


% 尝试从幅频图观察
fre = FRE_F(3,5);           % F调5
mpow = sortf(1:10,:);        % 取出各时刻最强的3个频率分量
[posR, posC] = find(abs(mpow-fre)<fre*0.05);    % 找到该音所在时域位置
if ~isempty(posC)
    
    fft_plot(fmt(T(posC(1))*Fs:T(posC(end))*Fs),'F调5');
    Fre_So = [1, 0.262, 0.156];
else
    Fre_So = 1;
end

fre = FRE_F(3,6);           % F调6
mpow = sortf(1:3,:);        % 取出各时刻最强的3个频率分量
[posR, posC] = find(abs(mpow-fre)<fre*0.05);    % 找到该音所在时域位置
if ~isempty(posC)
    
    fft_plot(fmt(T(posC(1))*Fs:T(posC(end))*Fs),'F调6');
    Fre_La = [1, 0.153, 0.168];
else
    Fre_La = 1;
end

fre = FRE_F(3,2);           % F调2
mpow = sortf(1:5,:);        % 取出各时刻最强的3个频率分量
[posR, posC] = find(abs(mpow-fre)<fre*0.05);    % 找到该音所在时域位置
if ~isempty(posC)
    
    fft_plot(fmt(T(posC(1))*Fs:T(posC(end))*Fs),'F调2');
    Fre_Re = [1, 0.423, 0.1194];
else
    Fre_Re = 1;
end

fre = FRE_F(3,1);           % F调1
mpow = sortf(1:5,:);        % 取出各时刻最强的3个频率分量
[posR, posC] = find(abs(mpow-fre)<fre*0.05);    % 找到该音所在时域位置
if ~isempty(posC)
    
    fft_plot(fmt(T(posC(1))*Fs:T(posC(end))*Fs),'F调1');
    Fre_Do = [1, 0.2621, 0.291];
else
    Fre_Do = 1;
end

fre = FRE_F(2,6);           % F调低音6
mpow = sortf(1:5,:);        % 取出各时刻最强的3个频率分量
[posR, posC] = find(abs(mpow-fre)<fre*0.05);    % 找到该音所在时域位置
if ~isempty(posC)
    
    fft_plot(fmt(T(posC(1))*Fs:T(posC(end))*Fs),'F调6-');
    Fre_La_ = [1, 0.9033, 0.216];
else
    Fre_La_ = 1;
end

dfh_score = [3, 5, 1;
    3, 5, 0.5;
    3, 6, 0.5;
    3, 2, 2;
    3, 1, 1;
    3, 1, 0.5;
    2, 6, 0.5;
    3, 2, 2];

mus = SoundSynthesis6_2_3_2(dfh_score, 120);
sound(mus, SampRate);

audiowrite("..\Output\ChineseRed6_2_3_2.wav", mus, SampRate);

% for i = length(funds)
%     pos = find(funds{i}==fre);
%     if ~isempty(pos)
%         break
%     end
% end
% 
% fprintf("频率为%f的基频大小为: %f\n",fre, fPs{i}(pos));

% figure
% basicFre = 329;
% index = find(f>basicFre*0.95&f<basicFre*1.05);
% plot(f(index), P1(index));
% [pks, locs] = enve(P1(index), 3);
% f1 = f(index);
% hold on
% xline(f1(locs));
% [M,I] = min(abs(f1(locs)-basicFre));

%% ConstantDefinition.m

function ConstantDefinition()
% ConstantDefinition 设置全局常量并保存在ConstantDefinition.mat里，用到时load即可

    % F调唱名对应频率
    % 行：大字组，小字组，小字一组，小字二组
    FRE_F = [
             87.31, 98, 110, 116.54, 130.81, 146.83, 164.81;
             174.61, 196, 220, 233.08, 261.63, 293.66, 329.63;
             349.23, 392, 440, 466.16, 523.25, 587.33, 659.25;
             698.46, 783.99, 880, 932.33, 1046.50, 1174.66, 1318.51;
             ];

    % used in function SingleSound
    PHI = 0;                % 初始相位
    TurnPoint1 = 0.25;      % 转折点
    TurnPoint2 = 0.35;
    TurnPoint3 = 0.7;
    Arate = 1.5;            % 最大值占平稳时幅值的比率

    HarmAmp1 = 0.2;         % 第一泛音列幅值占基音的幅值比例
    HarmAmp2 = 0.3;         % 第二泛音列幅值占基音的幅值比例

    % used in function SoundSynthesis
    SampRate = 44100;       % 采样率
    OverlapRate = 0.4;      % 重叠率，与前一个音有所重叠

    save("..\data\ConstantDefinition.mat", "FRE_F", "PHI", "TurnPoint1", "TurnPoint2", "TurnPoint3", "Arate", ...
        "SampRate", "OverlapRate", "HarmAmp1", "HarmAmp2");
end

%% 迭代计算包络
function [pks, locs] = enve(data, iter)
    if iter > 1
        [pks, l1] = findpeaks(data);
        [pks, l2] = enve(pks, iter-1);
        locs = l1(l2);
    else
        [pks, locs] = findpeaks(data);
    end
end

%% function SoundSynthesis
function res = SoundSynthesis6_2_3(musicScore, bpm)
% SoundSynthesis 将SingleSound得到的许多单音按一定的重叠比率整合为一条音频
%   musicScore: 乐谱参数，默认F调，以下方的fre_f为准。
%               size: [n 3] n行，1,2列分别为小字几组和对应唱名，3行为所占拍数
%               （后期可扩展，实现每个音单独设置其相应参数）（ps：什么UTAU）
%   OverlapRate: 重叠率，后一个音与前一个音的衰减部分比率
%   bpm: 每分钟拍数，衡量音乐速率
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[:,1]

    dur = 60 / bpm;     % 每拍所占时长（单位：秒）
    load("..\data\ConstantDefinition.mat", "FRE_F", "TurnPoint3", "SampRate", "OverlapRate");
    offset = 0;         % 重叠时向前偏移值
    offset_t = 0;       % 重叠时间（单位：秒）
    i = 1;              % 循环变量
    res = 0;
    for index = musicScore'
        sou = SingleSound6_2_3(FRE_F(index(1),index(2)), SampRate, 15, ...
                index(3) * dur + offset_t);
        i = i - offset;
        j = 1;
        while i <= size(res,1) && i > 1
            res(i, 1) = res(i) + sou(j);
            j = j + 1;
            i = i + 1;
        end
        while j <= size(sou,1)
            res(i, 1) = sou(j);
            j = j + 1;
            i = i + 1;
        end
        offset = ceil(OverlapRate * (1 - TurnPoint3) * size(sou, 1));
        offset_t = OverlapRate * (1 - TurnPoint3) * (index(3) * dur + offset_t);
    end
    res(i, 1) = 0;
end

%% fft_plot
function fft_plot(data, name)
    L = length(data);
    Y = fft(data);
    % 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
    P2 = abs(Y/L);
    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);
    
    f = 8000*(0:(L/2))/L;
    figure
    plot(f,P1) 
    title(strcat(name, " 频谱分析"));
    xlabel("f (Hz)");
    ylabel("|P1(f)|");
end

%% SingleSound
function res = SingleSound6_2_3(fre, SampRate, amplitude, durat, overtone)
% SingleSound: 单个乐音生成器（叠加泛音列）
% fre: 乐音频率
% SampRate: 采样率
% amplitude: 乐音幅值（持续时幅值）
% durat: 乐音时长（单位：秒）
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[时长*采样率 + 1， 1]

    res = SingleFre6_2_3(fre, SampRate, amplitude, durat);
    for i = 2 : length(overtone)
        res = res + overtone(i) * SingleFre6_2_3(fre*i, SampRate, amplitude, durat);
    end
end

function res = SingleFre6_2_3(fre, SampRate, amplitude, durat)
% SingleFre: 单频率乐音生成器
% fre: 乐音频率
% SampRate: 采样率
% amplitude: 乐音幅值（持续时幅值）
% durat: 乐音时长（单位：秒）
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[时长*采样率 + 1， 1]

    res = zeros(ceil(durat * SampRate), 1);
    load("..\data\ConstantDefinition.mat", "PHI", "TurnPoint1", "TurnPoint2", "TurnPoint3", "Arate");
    ind = 1;            % 循环变量
    for t = 0:1/SampRate:durat
        if t <= TurnPoint1 * durat
            % res(ind,1) = (Arate * amplitude / TurnPoint1 / durat * t) * sin(fre * t + PHI);
            res(ind,1) = (1-10^(-t))/(1-10^(-TurnPoint1 * durat)) * Arate * amplitude * sin(fre * t + PHI);
        elseif t <= TurnPoint2 * durat
            res(ind,1) = (amplitude*(Arate-1)*(exp(-(t-TurnPoint2*durat))-1)/(exp(-(TurnPoint1*durat-TurnPoint2*durat))-1)+amplitude)*...
                sin(fre * t + PHI);
            % res(ind,1) = (-(Arate - 1) * amplitude / (TurnPoint2 - TurnPoint1) / durat * t + ...
            %     (Arate * TurnPoint2 - TurnPoint1) * amplitude / (TurnPoint2 - TurnPoint1)) * ...
            %     sin(fre * t + PHI);
        elseif t <= TurnPoint3 * durat
            res(ind,1) = amplitude * sin(fre * t + PHI);
        else
            res(ind,1) = amplitude*(10.^((1/Arate)-(t-TurnPoint3 * durat)/(durat-TurnPoint3 * durat)*(1/Arate))-1)/(10.^(1/Arate)-1) * sin(fre * t + PHI);
            % res(ind,1) = (amplitude / (1 - TurnPoint3) / durat * (durat - t)) * ...
            %     sin(fre * t + PHI);
        end
        ind = ind + 1;
    end 
end

%% ffilter
function [fundamental, funP] = ffilter(f,fp)
% 在输入的一系列频率值中, 找出可能的基频
    fundamental = [];
    funP = [];
    for j = 1:length(f)
        if isempty(fundamental)
            fundamental = [fundamental,f(j)];
            funP = [funP,fp(j)];
        else
            ratio = f(j)./fundamental;
            inrange = (ratio<(round(ratio)*1.05)) + (ratio>(round(ratio)*0.95));    % 近似整数倍
            if all(inrange~=2)
                fundamental = [fundamental,f(j)];
                funP = [funP,fp(j)];
            end
        end
    end
end
 
%% SoundSynthesis.m
function res = SoundSynthesis6_2_3_2(musicScore, bpm)
% SoundSynthesis 将SingleSound得到的许多单音按一定的重叠比率整合为一条音频
%   musicScore: 乐谱参数，默认F调，以下方的fre_f为准。
%               size: [n 3] n行，1,2列分别为小字几组和对应唱名，3行为所占拍数
%               （后期可扩展，实现每个音单独设置其相应参数）（ps：什么UTAU）
%   OverlapRate: 重叠率，后一个音与前一个音的衰减部分比率
%   bpm: 每分钟拍数，衡量音乐速率
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[:,1]

    dur = 60 / bpm;     % 每拍所占时长（单位：秒）
    load("..\data\ConstantDefinition.mat", "FRE_F", "TurnPoint3", "SampRate", "OverlapRate");
    offset = 0;         % 重叠时向前偏移值
    offset_t = 0;       % 重叠时间（单位：秒）
    i = 1;              % 循环变量
    res = 0;
    for index = musicScore'
        sou = SingleSound6_2_3_2(FRE_F(index(1),index(2)), SampRate, 10, ...
                index(3) * dur + offset_t);
        i = i - offset;
        j = 1;
        while i <= size(res,1) && i > 1
            res(i, 1) = res(i) + sou(j);
            j = j + 1;
            i = i + 1;
        end
        while j <= size(sou,1)
            res(i, 1) = sou(j);
            j = j + 1;
            i = i + 1;
        end
        offset = ceil(OverlapRate * (1 - TurnPoint3) * size(sou, 1));
        offset_t = OverlapRate * (1 - TurnPoint3) * (index(3) * dur + offset_t);
    end
    res(i, 1) = 0;
end

%% SingleSound.m

function res = SingleSound6_2_3_2(fre, SampRate, amplitude, durat)
% SingleSound: 单个乐音生成器（叠加第一、二泛音列）
% fre: 乐音频率
% SampRate: 采样率
% amplitude: 乐音幅值（持续时幅值）
% durat: 乐音时长（单位：秒）
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[时长*采样率 + 1， 1]

    load("..\data\ConstantDefinition.mat");
    switch fre
        case FRE_F(3,5)
            harm = [1, 0.262, 0.156];
        case FRE_F(3,6)
            harm = [1, 0.153, 0.168];
        case FRE_F(3,2)
            harm = [1, 0.423, 0.1194];
        case FRE_F(3,1)
            harm = [1, 0.2621, 0.291];
        case FRE_F(2,6)
            harm = [1, 0.9033, 0.216];
        otherwise
            harm = 1;
    end
    res = 0;
    for i = 1:length(harm)
        res = res + harm(i) * SingleFre6_2_3_2(fre*i, SampRate, amplitude, durat);
    end
end

function res = SingleFre6_2_3_2(fre, SampRate, amplitude, durat)
% SingleFre: 单频率乐音生成器
% fre: 乐音频率
% SampRate: 采样率
% amplitude: 乐音幅值（持续时幅值）
% durat: 乐音时长（单位：秒）
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[时长*采样率 + 1， 1]

    res = zeros(ceil(durat * SampRate), 1);
    load("..\data\ConstantDefinition.mat", "PHI", "TurnPoint1", "TurnPoint2", "TurnPoint3", "Arate");
    ind = 1;            % 循环变量
    for t = 0:1/SampRate:durat
        if t <= TurnPoint1 * durat
            res(ind,1) = (Arate * amplitude / TurnPoint1 / durat * t) * sin(fre * t + PHI);
        elseif t <= TurnPoint2 * durat
            res(ind,1) = (-(Arate - 1) * amplitude / (TurnPoint2 - TurnPoint1) / durat * t + ...
                (Arate * TurnPoint2 - TurnPoint1) * amplitude / (TurnPoint2 - TurnPoint1)) * ...
                sin(fre * t + PHI);
        elseif t <= TurnPoint3 * durat
            res(ind,1) = amplitude * sin(fre * t + PHI);
        else
            res(ind,1) = (amplitude / (1 - TurnPoint3) / durat * (durat - t)) * ...
                sin(fre * t + PHI);
        end
        ind = ind + 1;
    end 
end